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An optimised algorithm for ionized impurity scattering in Monte Carlo simulations 

W.Th. Wenckebaclfl and P. Kinsleij] 

Department of Applied Physics, Technical University Delft, Lorentzweg 1, 2628 CJ DELFT, The Netherlands. 

We present a new optimised model of Brookes-Herring ionized impurity scattering for use in Monte Carlo 
simulations of semiconductors. When implemented, it greatly decreases the execution time needed for simula- 
tions (typically by a factor of the order of 100), and also properly incorporates the great proportion of small angle 
scatterings that are neglected in the standard algorithm. It achieves this performance by using an anisotropic 
choice of scattering angle which accurately mimics the true angular distribution of ionized impurity scattering. 



This paper was published as Computer Physics Communi- 
cations 143, 136 (2002). 



I. INTRODUCTION 

Ionized impurity scattering ||l|] dominates the transport 
properties of highly doped semiconductors, but accurate mod- 
eling of this scattering process is difficult because of the long 
range of the Coulomb force. The method in general use is 
that by Brookes-Herring, which is regarded as sufficiently ac- 
curate when simulating carrier mobilities in semiconductors. 
However, Brookes-Herring scattering is very anisotropic, so 
the standard overestimation-rejection algorithm[|| for imple- 
menting this in Monte Carlo simulations is hugely inefficient. 

The overestimation-rejection algorithm works as follows. 
The overestimated scattering rate is used to determine the free 
flight time of the particle, and at the end of every period of 
free flight a check is made to see if a scattering event has oc- 
curred. This, of course, leads to more checks being made than 
needed, but this problem is unavoidable since we cannot know 
in advance what the exact rate, for whatever final state that 
eventuates, will be. Each check involves three stages: selec- 
tion of a possible final state, a calculation of the scattering rate 
into that final state would be, and a decision as to whether the 
scattering occurs or not. 

This means that despite the fact that a scattering event does 
not occur at the end of every free flight, we still have to cal- 
culate the exact differential scattering rate for some particu- 
lar choice of scattering event. If the overestimation is too 
great, this means the program will calculate far more differ- 
ential scattering rates than strictly necessary; and then if this 
calculation is computationally expensive, the simulation will 
need to run for a much longer time. This is generally the case 
with ionized impurity scattering models in semiconductors, 
whose isotropic selection of possible final states means that 
overestimated rate must also be isotropic. Since the differ- 
ential scattering rate is strongly peaked for small angles, this 
means that the overestimated maximum scattering rate is far 
in excess of the typical differential rate. 
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Our new algorithm has the angular dependence of ionized 
impurity scattering built in, thus our overestimation becomes 
extremely efficient (see e.g. [H]). The resulting code therefore 
will run the same simulation faster than one with the stan- 
dard isotropic algorithm; or, alternatively, could run a much 
more accurate simulation in a comparable time. An early ver- 
sion of this algorithm was used recently to obtain manage- 
able execution times in simulations of hot-hole lasers in III-V 
semiconductors . 

The paper is organised as follows: section || describes the 
different algorithms; section |l| describes their comparative 
performance, and finally, in section [V we present our conclu- 
sions. 



II. IONIZED IMPURITY SCATTERING 

Here we are interested in optimising the algorithm for the 
Brookes-HeiTing model of ionized impurity scattering. This 
looks at the likelihood of the particle of interest (electron or 
hole) being deflected by the screened Coulomb field from a 
charged impurity [jl]]. The Brookes-Herring rate for ionized 
impurity scattering is 



W,iiff 



2% n,-,e^ 
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(1) 

Here n„ are the ionized impurity and hole (or electron) con- 
centrations, q-Q is the square of the minimum change of the 
fc-vector, ks is Boltzman's constant, T the temperature, EoE,- 
the permitivity, e the electron charge, / is the number of the 
initial band, / is the number of the final band. The so-called 
overlap factor is the square of the scalar product of the eigen- 
vectors corresponding to the initial and final state respectively, 

and is | (nfe; | mkf) |^. 

Note in particular that the rate W^/,// is proportional to 
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We would usually assume that go is just the inverse of the 
Debeye screening length, which is given by 
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although in subsection [I C we introduce an additional con- 
tribution k,„i„, where ^^q^+kl^i„- 

The difficulty with ionized impurity scattering is that the 
scattering rate becomes very large for small values of ^ = 
|k/ — k, |. We can distinguish two situations leading scattering 
rates that are very large compared to the average scattering 
rate: 



(1) small angle scattering, i.e. q is small compared to the 

absolute value of the incoming wave vector kj = |k, |. 
Note that ki is of the same order of magnitude as the 
absolute value of the outgoing wave vector A;/ = |k/|. 

(2) low hole or electron energy, i.e. where ki is small by 

itself. 
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Qi = 2sin(e/2) = ^2(1 -cose) 



(4) 



(5) 



where 9 is the scattering angle. We now introduce an 
anisotropic weighing function /(9,(|)) for the probability of 
finding the direction of the final k-vector. This same factor 
will then be used to compensate for this weighting by subse- 
quently dividing the differential scattering by /(9,(|)). At first 
sight one would be inclined to choose 



A. Standard Isotropic Algorithm 

The great variation in the scattering rate by angle and wave 
vector means that it is not straightforward to efficiently model 
the scattering process. If we take the simple approach of an 
isotropic choice of scattering angle, the usual algorithm over- 
estimates the differential scattering rate W^'^^ 9, (|)) by 
its maximum value W,'u[/ . Therefore our Monte Carlo pro- 
gram will select an ionized impurity scattering for any parti- 
cle at the overestimated rate WmH , then reject those for which 
the combination of exact scattering rate Wp''^'^(£'<^(f),9,(|)) and 
choice of random number do not lead to a scattering. 

To choose a possible final state for an ionized impurity scat- 
tering process, we start by choosing a random direction for 
k/ with all possible final directions have the same statistical 
weight /(9,(|)) = l/(47i). Once we have chosen this direc- 
tion, we can then work out the final k vector, and hence the 
true scattering rate. We compare this true rate to the overesti- 
mated differential scattering rate, and use the normal rejection 
procedure to determine statistically whether the scattering has 
occurred or not. 

For ionized impurity scattering Wma is much and much 
larger than Wp''^^(£'<^(f),9,(|)) for about any angle; and also 
this maximum value Wm^/ depends sensitively on the screen- 
ing length. This means that Wmlx^ is an inefficient overesti- 
mation, and as a result we do many calculations which do not 
lead to a scattering process. This significantly lengthens the 
execution times for our simulations. 



B. New Anisotropic Algoritlim 

We have improved significantly on the execution times of 
simulations by introducing an anisotropic choice of the scat- 
tering angle, whilst compensating for this by adjusting the fi- 
nal scattering rate. The angular dependence of the scattering 
rate is centered on the direction of the initial k vector, so it is 
useful to transform to angular coordinates (9,(|)) centered on 
this direction. For this purpose we first define 



Q\ 16sin4(9/2) 4(l-cos9)2 



(6) 



but this function cannot be normalised properly because the 
integral of /(9, (|)) over 9 and (|) diverges. We could try to solve 
this problem by assuming a minimum scattering angle 9o, so 
we get a minimum value for Qi equal to Qq = 2sin(9o/2) 
and to choose for Qq e.g. the approximate Debeye-value of 
Ridley[|lll, we have 



(60^,) 



2 



(7) 



We would then assume 9o to be so small that we may as- 
sume the absolute value of the incoming and outgoing wave 
vectors to be the same. However, excluding scatterings with 
q < qo leads to erroneous results because it is exactly these 
small angle scatterings which dominate ionized impurity scat- 
tering. Test simulations confirm that the normalisation proce- 
dure has an equally important influence on the final result as 
the way screening is implemented. 

To avoid erroneous results we choose /(9,(|)) to have the 
Brookes-Herring shape (see eqn. ||) - 
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{Q\ + Qlf (4sin2(9/2) + e2)' (2(l-cos9) + e2) 



The function can be normalised because 
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is finite. As a result the normalised statistical weight of 
finding a given direction is given by 



Qoi^+Qo) 



Go(4 + Gg) 



(Qi + Ql) 47t[2(l-cos9) + G2]^ 

(10) 

The next problem is to enter this statistical weight in the 
algorithm. To determine 9, (|) we choose two random numbers. 
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ri and r2, between and 1. In the case of an isotropic choice 
of direction we would have /(9,(|)) — l/4-n, and we would 
choose 



(11) 

and r2 = 2%{l -cosB) f(e,<^). (12) 
For our new algorithm, we choose ri as in eqn. (|ll]), and 

r2 = 7t[2(l-cose) + e2]/(e,^)_M0 (13) 

The extra terms with Qq are needed to ensure that r2 lies 
between and 1 , so that — 1 < cos 9 < 1 . The weighed random 
choice of (|) and 9 is now obtained by solving — 2k ri and 



Note that this is the point where choosing a good way to 
normalise the weighing function /(9,(|)) pays off. The result- 
ing function becomes extremely simple! Finally, the differen- 
tial density of states is replaced by the total density of states 
divided by 4jt. Thus the overestimated scattering rate is given 
by 



W,o, = — 



2n nue'^ 



(19) 



Almost always ^ q^. 



so this overestimation is orders of 
magnitude smaller than that from the standard isotropic al- 
gorithm, leading to far fewer calculations of the differential 
scattering rate that are followed by rejection. 



C. Modified Anisotropic Algoritiim 



cos 9 



l.f 



Ql 



(14) 



Note that cos 9^1 for r2 ^ 1 , so small angle scattering is 
not excluded. 

In the new algorithm we need to compensate for the higher 
occurrence of small scattering angles caused by the introduc- 
tion of a weighing function /(9,(|)) ^ 1 /47t by multiplying the 
differential scattering rate by 



471 
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We construct an overestimated differential scattering rate 
by starting to look for the largest differential ionized impurity 
scattering rate including the factor /(9,(|)). The overlap factor 
is always less than one. The term 



1 



can be overestimated by 

C 



(16) 



{Ql 



(17) 



where C is equal to 1 if the band is isotropic, but may need 
to be larger for anisotropic bands. In practice, C needs to be 
determined by computer experiment, which in our simulations 
was set equal to 1.2. There is also our new weighting factor 
/(9,(|)), which with the previous factor (eqn. |l5|) may both be 
overestimated by 



Ql) 



c 



c 



c 



eg (4+ eg) (e?+eg)2^f 



eg (4 + eg) ql{4kj+ql) 
(18) 



Ionized impurity scattering is an elastic process, which 
means that small changes of the ^-vector imply a small change 
of the state of the hole or electron. A minimum change of 
the A;-vector follows immediately from the implementation of 
the Debeye screening through eqns. ([l|) and (|^). The result- 
ing minimum value of qo is proportional to ^/NH, so it be- 
comes very small for low ionized impurity concentrations. As 
a result, despite the anisotropic scattering angle selection, the 
fraction of overestimations and hence the execution times of 
the simulation will still become very large for low concentra- 
tions. We can therefore improve the overestimation by intro- 
ducing a second criterion determining whether the change of 
the fc-vector is small or not. Now the average energy of an 
electron or hole is given by 



3 h^<k^> 

— kgT = 

2 2m<,m* 



so the average value of k is given by 



(20) 
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k'> = 



3mem*kBT 



(21) 



This induces us to propose that we consider k small if its 
square is much smaller than this value. So we set 



D 



3t7ietn*kBT 



(22) 



where we choose D to be a small number. In practice we 
find that D = 0.01 is acceptable. We can now set the minimum 
change in wave vector to be 



1Q 



Id 



(23) 



Note that excluding scatterings with a 'small' change Ak 
of k means skipping an infinite number of small angle scat- 
terings. However these excluded scatterings have a negligible 
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influence on the mobility since the motion of the carrier is 
largely unaffected. E.g. if < cko, where k^/2m* — kyT 
and c = 0.1, then on the average such a scattering changes the 
component of k in the direction of the electric field by less 
than 1 %. To verify this argument we repeated a number of the 
simulations as varying c from 0.1 to 1.6. For c < 0.3 we saw 
no significant effect on the calculated mobilities. 



III. RESULTS 

To check the relative efficiency of these algorithms, we 
took an existing Monte Carlo code (used in Dijkstra and 
Wenckebach!^) ™d rewrote one of the subroutines in two 
ways: first, in accordance with the isotropic algorithm, and 
second, according to the modified anisotropic algorithm. We 
could then compare the performance of the two algorithms 
by comparing the performance of the two codes compiled 
with either of the two coded algorithms. In all other aspects 
the two codes were the same. However, with the isotropic 
choice of the scattering angle, small angle scatterings such 
that sin 9 < had to be excluded in order to let the program 
run within finite time. 

Both the material parameters we used, and the general 
method followed by the code were the same as those of Hinck- 
ely and Singh The codes were compiled with GNU 

g77 and run on a 500 MHz Pentium III running RedHat Linux 
6.2. The calculation was of the drift velocities of holes in un- 
strained silicon, where the magnetic field was equal to zero, 
the electric field equal to 5 kV/cm and the temperature 300 K. 
Only the fc-vector was integrated during free flight. Each sim- 
ulation consisted of either 25 blocks of 1000 real scatterings 
each or of 2 blocks of 100 real scatterings. Only in the former 
set of simulations could mobilities be determined. 



TABLE I: A comparison of execution times and mobilities for the 
isotropic and anisotropic algorithms. The material system was un- 
strained silicon with and electric field equal to 5 kV/cm and a temper- 
ature of 300 K. Note that in almost all cases, the simulations averaged 
over 25 x 1000 scatterings take hours using the isotropic algorithm, 
but only minutes using the anisotropic one. 



Code: 


isotropic choice of 


anisotropic choice of 


Scatters. 


25 X 


1000 


2 X 100 


25 X 


1000 


Cone. 


Exec, time 


Mobility 


Exec, time 


Exec, time 


Mobility 


cm^-' 


hh.mm.ss 


cm^A's 


hh.mm.ss 


hh.mm.ss 


cm^/Vs 


10'^ 


01.18.35 


167 ± 26 


00.00.38 


00.02.47 


95 ±40 


1018 


13.23.44 


210 ±46 


00.07.39 


00.02.17 


198 ± 50 


1017 


± 145 h 




01.09.47 


00.02.16 


386 ± 75 


loif- 


± 380h 




03.04.02 


00.03.40 


400 ± 60 


1015 


± 80h 




00.38.43 


00.06.39 


412 ± 52 


1014 


09.56.23 


452 ± 32 


00.05.07 


00.07.41 


430 ± 43 


1013 


01.06.30 


418 ±38 


00.00.32 


00.07.55 


420 ± 37 





00.06.37 


408 ± 50 


00.00.03 


00.06.37 


408 ± 50 



In Table | the first column gives the concentration of ionized 
impurities, the next three columns are data from the isotropic 
code, and the last three columns are data from the anisotropic 
code. The errors indicate the 95 % certainty level and corre- 
spond to about twice the standard deviation. In the case of an 
isotropic choice of the angle after scattering, with 25 blocks 
of 1000 real scatterings execution time became too long for 



10'^ to 10" cm 



Therefore it was estimated from 



runs with 2 blocks of 100 scatterings. We can see that in all 
cases the anisotropic method was significantly faster than the 
isotropic - the longest anisotropic execution time is 7 minutes 
55 seconds, whereas all but the n,-, = simulation using the 
isotropic code took over 1 hour. 



TABLE II: A comparison of execution times and mobilities for sim- 
ulations with the new anisotropic algorithm, with comparison to the 
tabulated mobilities from Sze[8]. These are for unstrained silicon at 
a temperature of 300 K and are averaged over 50 x 5000 scatterings. 



Cone. 


Exec, time 


Mobility (simul.) 


Mobility (Sze) 


-3 

cm 


hh.mm.ss 


cm^A^s 


cm^A's 


1015 


00.27.50 


132 ± 13 




1018 


00.22.18 


196 ± 21 


150 


1017 


00.22.34 


334 ± 24 


300 


IQl^ 


00.36.29 


400 ± 19 


430 


1015 


01.06.43 


418 ± 14 


460 


1014 


01.17.14 


422 ± 14 


470 


1013 


01.19.12 


430 ± 12 







01.02.28 


428 ± 14 





As a final check, table |l| shows simulation results compared 
with experimental values for the mobility in silicon as found 
in Sze[^. For more precision each simulation consisted of 50 
blocks of 5000 real scatterings, and the electric field was kept 
at 5 kV/cm. Note that again the stated errors correspond to 
a 95 % certainty level, and that we have good agreement be- 
tween our code with the optimised ionized impurity scattering 
algorithm and the values given by Sze. 



IV. CONCLUSIONS 

We have demonstrated a new way of significantly reducing 
execution times in simulations of systems involving Brooks- 
Herring ionized impurity scattering. This was done by in- 
troducing a carefully calculated anisotropic scattering angle 
for these ionized impurity processes, thus avoiding the usual 
problem of a large fraction of inefficient overestimations. Ex- 
ecution times of hole mobility calculations in silicon showed 
that the speedup was strongly dependent on impurity concen- 
tration, and generally well in excess of a factor of 5, typically 
being of the order of 100. The method is applicable to any 
simulation in which Brookes-Herring ionized impurity scat- 
tering is implemented, and it could be usefully generalised to 
other strongly anisotropic scattering processes, such as polar 
optical phonon and piezoelectric scatterings. 
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